A fundamental problem with MRI in prior art is that RF emitting nuclei are distributed in a 3D volume space (x,y,z) but only one fully parallel spatial encoding scheme is available, that is frequency encoding. The phase encoding scheme is only partly parallel. Therefore sequential scanning is needed along at least one of the three dimensions (x,y,z) to image a target object in 3D space. This makes MRI in prior art slow. A fully parallel MRI to image a 3D volume space would need 3 independent fully parallel encoding schemes, one for each of the three Cartesian coordinates (x,y,z). In principle, the present invention FIT provides this facility. The three-dimensional spatial distribution f(x,y,z) of an RF emission source is encoded by the coordinates (x,y,z) themselves, in the sense that this f(x,y,z) uniquely determines and vice versa the 3D intensity distribution field g(x,y,z) of the RF emission produced by it in its vicinity in a 3D volume space around it. In other words, the 3D intensity distribution field g(x,y,z) is a function of (or depends on) three independent parameters, i.e. x,y, and z, not just on only one or two parameters such as frequency and phase in conventional MRI. FIT can be used as a stand-alone scheme for spatial encoding, or it can be combined with frequency encoding scheme for parallel MR imaging. It can also be combined with phase encoding scheme, or both frequency encoding and phase encoding schemes together. This technique of combining FIT with frequency and/or phase encoding schemes permits different trade-offs in the 3D imaging of objects in terms of speed, accuracy, cost, machine size, etc. In the preferred embodiment two-dimensional (2D) planar slices of a 3D target object are scanned sequentially along the Z-axis, and within a 2D slice location along the X-axis is encoded by frequency and Y-position is encoded by the Y-coordinate itself or the dependence of intensity or amplitude of radiation on the Y-coordinate. In this case, the conventional phase encoding has been replaced by spatial-coordinate encoding or amplitude encoding. Unlike the conventional phase-encoding technique, the method of the present invention is completely parallel and therefore much faster.
Principles of Field Image Tomography (FIT) in a different form with no name was first developed and applied for SPECT and PET imaging in the following patent application by this inventor, the teachings of which are incorporated herein by reference:    M. Subbarao, “Method and Apparatus for High-Sensitivity Single-Photon Emission Computed Tomography”, U.S. patent application Ser. No. 12/586,863, filed on Sep. 29, 2009.
In the present invention FIT approach has been modified in a novel way and applied to Magnetic Resonance Imaging. A good introduction to an exemplary MRI system and associated problems can be found in the following US patent (U.S. Pat. No. 6,476,606 B2), the teachings of which are incorporated herein by reference:    R. F. Lee, “Method for Parallel Spatial Encoded MRI and Apparatus, Systems and Other Methods Related Thereto”, U.S. Pat. No. 6,476,606 B2, Nov. 5, 2002.
In particular, MRI systems in prior art are very slow due to the so called phase-encoding step in data capture. This phase-encoding step can be avoided in the present invention for fast data capture and image reconstruction. As mentioned earlier, FIT can be used on its own or it can be combined with frequency or/and phase encoding techniques to achieve different trade-offs in speed, accuracy, cost, etc.
The method and apparatus of the present invention are specifically related to exploiting the information in the spatial intensity distribution of the three-dimensional (3D) radio frequency (RF) emission field in a 3D volume space. This information can reveal the 3D spatial density distribution of the RF emission source. In particular, the 3D volume space, when necessary, is designed to extend substantially along the radial distance pointing away from the RF emission source. The intensity of RF emission decreases with the square of the radial distance from the source and this variation in 3D space provides information on the location of the source. Therefore measuring the RF emission along different distances along the radial direction in addition to different angular directions is important.
FIT is based on analyzing the emission field intensity distribution in a three-dimensional (3D) volume space around the emission source. In particular, radiation intensity from a point source decreases with increasing square of the distance from the source, and this characteristic of radiation propagation is exploited in determining the location and intensity of point sources by measuring radiation intensity at different radial distances and analyzing the data. The radial distance r depends on all the three spatial coordinates x,y, and z, and therefore spatial localization in all the three coordinates becomes possible.
MRI machines in prior art are slow due to sequential scanning of different slices of the target object at different time instants. They do not include modeling and exploiting or measuring RF emission in a 3D volume space that extends substantially along the radial direction, even when that is needed and would be helpful.
Although there are a variety of spatial encoding methodologies or techniques being implemented, the most popular method used in commercial MRI scanners is two dimensional Fourier transform (2DFT) encoding in which a two-dimensional spatial plane (e.g., XY plane) is encoded with both frequency and phase of the MR signals. Typically during one data acquisition, only a one dimensional time-domain signal is obtained and thus 2DFT encoding requires repeating the data acquisitions sequentially to achieve a pseudo second dimension of the time domain signals. The second dimension of the spatial information is encoded into the phase component by repeating the data acquisition with different phase encoding gradient strengths i.e., varying magnetic field gradient along the y-axis Gy to create the other pseudo-time dimension. In this technique, each slice of the imaged object is in effect divided into a multiplicity of layers in the y-direction or axis corresponding to the number of pixels in that direction (e.g., 128, or 256). The number of pixels in turn is representative of the desired image resolution. The x-direction scanning process or the data acquisition is repeated by sequentially and repeatedly stepping through each of these y-direction layers. Because the resolution of the time-domain signals depends on the number of repetitions of the data acquisitions, and the repetition rate is limited by the MR relaxation times, a higher resolution image takes a longer time. In the present invention, this limitation can be relaxed by using more number of Radio Frequency Sensor Coils (RFSCs) at different positions and distances, and therefore the present invention can be made much faster.
There is a tradeoff between spatial resolution and imaging time because higher resolution images require a longer imaging time. This balance between spatial and temporal resolution is particularly important in cardiac MR, where fine details of coronary artery anatomy must be discerned on the surface of a rapidly beating heart. Thus, a high-resolution image acquired over a large fraction of the cardiac cycle, will be blurred and distorted by the motion of the beating heart.
In one embodiment of the present invention, the phase encoding is avoided by measuring the RF emission from one column of voxels parallel to the y-axis simultaneously with a plurality of RF sensitive coils (RFSCs) that are arranged in a 3D volume space in a non-degenerate configuration or arrangement. All voxels in a given column emit RF radiation at a constant frequency, but voxels in different columns emit RF radiation at different frequencies. The time-dependent emission radiation is measured by each coil and the resulting one-dimensional time signal is Fourier transformed to obtain the magnitude of different frequencies. The magnitude at each frequency gives the total radiation intensity emitted by all voxels in the corresponding column. The frequency itself encodes the position of the corresponding column along the x-axis. Different RFSCs are at different positions and distances from voxels in a column and therefore record different magnitudes for a given frequency. This measured difference in the magnitude of frequency components ultimately reveals the structure of the distribution of emission intensity from different voxels in a given column. The number of RFSCs must be comparable to the number of voxels in a column. They may lie on a line, curve, planar surface, curved surface, or distributed in a 3D volume space when necessary, so that the configuration is non-degenerate, i.e. the configuration is such that it captures and provides all the available and necessary information.
The present invention paves the road for true parallel MRI and the achievement of manifold reductions in minimum MRI scan-time for rapid dynamic studies that require high time resolution. The present invention can be used to find the 3D tomography of any field emitting source. For example, as mass emits gravitational field, the 3D density of earth, moon, and other such bodies can be determined by measuring the gravitational field in a large 3D volume space around these objects. The 3D volume space, when necessary, is extended both in radial and angular directions. The basic principle remains the same as that described in the present invention. The present invention can also be used with radiation frequencies (e.g. millimeter waves or microwaves) other than Radio Frequencies along with suitable magnetic fields. The present invention can also be used to make small (about the size of a hair dryer or a kitchen microwave) portable FIT machines for MRI.
2.1 Theoretical Basis of the Present Invention
The insight and inspiration that led to the present invention are two inventions for image deblurring and 3D shape-from-defocus techniques disclosed by this inventor in the following two patents:    1. M. Subbarao, “Methods and apparatus for computing the input and output signals of a linear shift-variant system”, Jul. 7, 2009, U.S. Pat. No. 7,558,709.    2. M. Subbarao, “Direct Vision Sensor for 3D Computer Vision, Digital Imaging, and Digital Video”, Aug. 18, 2009, U.S. Pat. No. 7,577,309.
The present invention is based on a new theory. It is based on measuring and inverting weighted volume integrals instead of simple line integrals corresponding to Radon Transform used in the past in older MRI methods. The weighting factor in the volume integrals corresponds to the reciprocal of the square of the distance between a field emitting source element and the field measuring element (RFSCs in MRI). In MRI, RF sensitive coils are used to measure total emission from many or all the emission source elements in a certain 3D volume space.
Another fundamental feature of the new theory is that weighted Volume integrals must be measured in a 3D volume space that may extend substantially along the radial direction pointing away from the emission source. Measuring on a thin surface that is roughly perpendicular to the direction of emission rays as in the conventional MRI theory may not be adequate in some problems. For example, RF emission from a source inside a human body may be measured in a 3D volume space in the shape of a thick annular cylinder, with an inner radius of 200 mm and an outer radius of 600 mm. Alternatively, measurements can also be taken on a small set of concentric circles of different radii, or concentric polygons of different sizes. Therefore, in the present invention, weighted 3D volume integrals of spatial density distribution of the emission source are measured at a set of points that could be spread out in a 3D volume space.
The validity of the present invention has been verified through simple computer simulation experiments. A discrete source distribution was generated using a random number generator. The measured radiation field due to this distribution was computed at a discrete set of points by using a model of radiation propagation. This measured data was successfully inverted to obtain the original spatial distribution of the radiation source.
2.2 Detailed Theory: Deriving H and g=Hf+n
In this section additional details are provided on deriving a system matrix H and the equation g=Hf+n which is the basis of the method of the present invention.
A 3D volume space V is selected in a target object whose 3D tomography needs to be imaged by MRI. This volume V of the target object is placed in a known magnetic field B. In the first embodiment of the present invention, B is a constant field within volume V. The nuclei of atoms in the material substance of the target object are usually spinning like a top and produce a magnetic field of their own. These spinning nuclei will become aligned so that their magnetic fields will be parallel or anti-parallel with the magnetic field B. In addition, magnetic field B exerts a torque on the spinning nuclei, and therefore they precess like a spinning top with Larmor frequency around the axis along the direction of B. Larmor frequency w is given by w=γB where γ is the gyromagnetic constant of the nuclei which is different for nuclei of different elements.
A suitable Radio Frequency (RF) pulse is beamed into the volume V of the target object so that many of the nuclei whose magnetic field are parallel or aligned with B are excited through energy absorption and become anti-parallel to B. (A second RF beam may be sent to realign the magnetic axes of the nuclei for some purpose such as facilitating measurement.) These excited nuclei then gradually lose their absorbed energy and return to their normal state which is parallel to B, and in this process the lost energy is emitted as RF emission. The amount or intensity of RF emission from each small volume element or voxel in V is a characteristic of the density of different types of nuclei, e.g. Hydrogen nuclei of water, in the material of the target object in V. By determining the intensity f(x,y,z) of this RF emission from each voxel located at Cartesian coordinates (x,y,z) in volume V of the target object, the 3D tomographic Magnetic Resonance Image f of the target object is obtained. However, in practice, it is possible to measure only the total sum or volume integral of RF emission from all the voxels in V. It is not possible to measure the RF emission intensity f(x,y,z) separately from individual voxels at (x,y,z). Therefore this measured volume integral (which is weighted by the factor of reduction in radiation intensity with distance from the source element) needs to be processed to reconstruct the original image f(x,y,z).
Radio Frequency Sensitive Coils (RFSCs) are used to measure the sum or volume integral of RF emission from all the voxels in V. In the prior art on MRI, this problem of measuring emission from individual voxels is solved in a few different ways. In the first way, a method similar to that in X-ray Computed Tomography is used. In this method, volume V is restricted to be a thin line or column of volume space. This is achieved by using a magnetic field B that linearly varies along one dimension, say the Z-axis in a Cartesian coordinate system centered in the target object. Then an RF pulse of a certain frequency is sent to select and excite only a thin two-dimensional (2D) slice or cross-section parallel to the X-Y plane in the target object that is perpendicular to the Z-axis. Within this 2D plane, a single line or column of voxels is frequency encoded by introducing a gradient magnetic field perpendicular to the Z-axis, say a field that increases linearly along the positive X-axis. In this case the field is a constant along lines parallel to the Y-axis and therefore voxel columns parallel to the Y-axis contain nuclei that are spinning with a constant frequency. On the other hand, the frequency of nuclei will increase from one column to the next along the positive X-axis. In this way spatial location along the X-axis of columns of voxels are frequency encoded by associating a unique frequency for each column position along the X-axis.
RF coils sense the total emission signal from all the voxels in the selected slice that has been excited by the RF pulse, and the Fourier Transform of the sensed signal gives a parallel projection or line integral of the RF emission from the slice. By rotating the direction of the gradient field around the Z-axis (to say 90 different angles at 2 degree intervals in the range 0 to 180 degrees) such a projection (also called the Radon Transform) of the slice can be obtained for (around 90) different angles in the X-Y plane. Then, a tomographic image reconstruction algorithm such as the Filtered-Back Projection (FBP) or Algebraic Reconstruction Technique (ART) can be used to determine the two-dimensional (2D) intensity distribution f(x,y) of the RF emission from voxels in the selected slice of the target object. This method is repeated for different thin successive 2D slices perpendicular to the Z-axis, and the results are stacked together in proper order to obtain a 3D tomographic image of the entire target object.
The above method is slow due to the scanning of different 2D slices in sequence at first, and then scanning of different angles in each slice. This scanning may take up to 30 minutes for a whole body scan of a target object during which a patient may have to be sedated to avoid motion of the body.
In another approach, instead of rotating the gradient magnetic field to different angles in the previous method, spatial position of voxels within a column are further encoded by a complicated phase encoding technique, and a 2D Fourier Transform is used for tomographic image reconstruction. This phase encoding technique involves sending RF pulses repeatedly and taking measurements. This method is faster than the first approach, but still much slower than the present invention. The present invention is estimated to take around 80% less time than this phase-encoding method as there is no sequential scanning. Instead, a plurality of RFSCs placed at different positions and distances are used to simultaneously measure the radiation field intensity in a 3D volume space around the target object. This measured data is then processed to reconstruct the 3D tomography of the target object.
In the present invention, in principle, there is no need to scan either one 2D slice at a time in a 3D volume, or scan projections at different angles or phase-encode within a 2D slice. In the present invention, 3D tomography can be obtained at once by taking measurements of RF emission from the entire 3D volume V at once. However, in practice, scanning one 2D slice at a time, but reconstructing an entire 2D slice at once without scanning for angles/directions or phase is a good approach in the method of the present invention. Otherwise the number of RFSCs needed will be too many and this poses difficulties in designing the imaging machine. It is also possible to use both frequency encoding and phase encoding to select a volume space V, and then within this volume space apply the method of the present invention to determine higher resolution tomographic image of the target object.
Next the equation g=Hf+n will be derived for the field image g. An RF emission source Q of unit intensity (power or energy per unit time, e.g. 1 watt) at a point (x,y,z) as in FIG. 4 produces an emission field of intensity 1/(4πr2) (watts/meter sq.) at a distance of radius r from the source. This is due to the fact that the emitted 1 unit of power is uniformly distributed over a spherical surface of area 4πr2 as the energy propagates from the source point to the surface of the sphere.
Let f(x,y,z) be the density (watts/meter cubed) of the RF source at a point Q with coordinates (x,y,z) in a body. The RF source element in a small volume element dv=dxdydz at Q(x,y,z) produces an RF emission field at another point P(x′,y′,z′) as shown in FIG. 4 given byf(x,y,z)dxdydz h(x′−x,y′−y,z′−z) watts/m^2  (Eq. 2.1)whereh(x′−x,y′−y,z′−z)=1/(4λr2)  (Eq. 2.2)andr2=(x′−x)2+(y′−y)2+(z′−z)2  (Eq. 2.3)
Let a small sensor element of area dA with unit surface normal vector k pointing away from the RF source be placed at (x′,y′,z′). Let the angle between this vector k and direction of RF emissioni=(x′−x,y′−y,z′−z)/|(x′−x,y′−y,z′−z)|be θ(x′,y′,z′,x,y,z). For example, if area dA is facing the origin withk=(x′,y′,z′)/|(x′,y′,z′)|then
                              Cos          ⁡                      (                          θ              ⁡                              (                                                      x                    ′                                    ,                                                            y                      ⁢                                                                                                            ′                                    ,                                      z                    ′                                    ,                  x                  ,                  y                  ,                  z                                )                                      )                          =                ⁢                  k          ·          i                                        =                ⁢                                            x              ′                        ⁡                          (                                                x                  ′                                -                x                            )                                +                                    y              ′                        ⁡                          (                                                y                  ′                                -                y                            )                                +                                                    z                ′                            ⁡                              (                                                      z                    ′                                    -                  z                                )                                      /                                                          ⁢                  (                                                                  (                                                                            x                      ′                                        -                    x                                    ,                                                            y                      ′                                        -                    y                                    ,                                                            z                      ′                                        -                    z                                                  )                                                    ⁢                                                        (                                                      x                    ′                                    ,                                                            y                      ⁢                                                                                                            ′                                    ,                                      z                    ′                                                  )                                                            )                    
For a given RF emission frequency u, the RF emission intensity at the 3D point P(x′,y′,z′) due to all of the source f(x,y,z) emitting RF at frequency u is the sum or integral of the field due to each of the volume element, given by:
                              g          ⁡                      (                                          x                ′                            ,                                                y                  ⁢                                                                                        ′                            ,                              z                ′                                      )                          =                                            ∫                              ∫                ∫                                      V                    ⁢                      h            ⁡                          (                                                                    x                    ′                                    -                  x                                ,                                                      y                    ′                                    -                  y                                ,                                                      z                    ′                                    -                  z                                            )                                ⁢                      f            ⁡                          (                              x                ,                y                ,                z                            )                                ⁢                                          ⁢                      Cos            ⁡                          (                              θ                ⁡                                  (                                                            x                      ′                                        ,                                                                  y                        ⁢                                                                                                                      ′                                        ,                                          z                      ′                                        ,                    x                    ,                    y                    ,                    z                                    )                                            )                                ⁢                      ⅆ            x                    ⁢                      ⅆ            y                    ⁢                      ⅆ            z                                              (                  Eq          .                                          ⁢          2.4                )            The above equation is named the Field Image Equation (FIE). In the above equation, only the intensity or amplitude of the RF emission at a single frequency u is considered. However the time dependent signal recorded by RFSCs may contain a plurality of frequencies if FIT is combined with frequency encoding for one of the dimensions, say along the x-axis. In this case the time signals are Fourier transformed and the magnitude of the Fourier transform at frequency u is then used to obtain one equation as above. For each different frequency u, one separate equation as above is obtained and solved. Therefore g(x′,y′,z′) in the above equation should be considered as the magnitude of the one-dimensional Fourier Transform of the actual time-dependent signal g′(x′,y′,z′,t) at some frequency u wherein the Fourier Transform is computed with respect to the time variable t. In this case, a one-to-one correspondence exists between the frequency u in the Fourier domain and the spatial position x of f(x,y,z) which has been frequency encoded by u. Therefore we use g(x′,y′,z′) to denote the measured amplitude at (x′,y′,z′) for a given frequency u corresponding to x in f(x,y,z).
If, as in FIG. 4, a small sensor element P is placed at (x′,y′,z′) with unit surface area dA and the surface normal is at an angle of θ(x′,y′,z′,x,y,z) with respect to the direction of incidence of RF emission from the source volume element at (x,y,z), then the measured RF field intensity g(x′,y′,z′) at (x′,y′,z′) isg(x′,y′,z′)=f(x,y,z)h(x′−x,y′−y,z′−z) dxdydz cos(θ(x′,y′,z′,x,y,z)) watts/m^2  (Eq. 2.4).(but the actual measured power by one sensor element is g(x′,y′,z′) dA watts).In MRI, attenuation of RF emission as it traverses through the target object is taken to be negligible.
If this emission is attenuated by the RF sensitivity coil geometry by a factor of a(x′,y′,z′,x,y,z), then the measured radiation will beg(x′,y′,z′)=f(x,y,z)h(x′−x,y′−y,z′−z) dxdydz cos(θ(x′,y′,z′,x,y,z)) a(x′,y′,z′,x,y,z)  (Eq. 2.6)If there is no such attenuation, then a(x′,y′,z′,x,y,z)=1 for all (x′,y′,z′,x,y,z).
In the discrete domain let us define the system matrix H(x′,y′,z′,x,y,z) asH(x′,y′,z′,x,y,z)=h(x′−x,y′−y,z′−z)cos(θ(x′,y′,z′,x,y,z))a(x′,y′,z′,x,y,z)  (Eq. 2.7).Therefore, the measured RF field intensity (amplitude) due to the presence of all RF source elements distributed in the (x,y,z) space with density f(x,y,z) can be expressed in the discrete domain as
                              g          ⁢                      (                                          x                ′                            ,                                                y                  ⁢                                                                                        ′                            ,                              z                ′                                      )                          =                              ∑            x                    ⁢                                    ∑              y                        ⁢                                          ∑                z                            ⁢                                                          ⁢                                                H                  ⁡                                      (                                                                  x                        ′                                            ,                                                                        y                          ⁢                                                                                                                                ′                                            ,                                              z                        ′                                            ,                      x                      ,                      y                      ,                      z                                        )                                                  ⁢                                                      f                    ⁡                                          (                                              x                        ,                        y                        ,                        z                                            )                                                        .                                                                                        (                  Eq          .                                          ⁢          2.8                )            In the above equation, unit volume was taken to be that of one voxel so that the term dxdydz=1. The summation is carried-out over all volume elements or voxels at points (x,y,z) where the RF source may be present. Due to noise and measurement errors together contributing n(x′,y′,z′) to the above measured value, we obtain:
                              g          ⁢                      (                                          x                ′                            ,                                                y                  ⁢                                                                                        ′                            ,                              z                ′                                      )                          =                                            ∑              x                        ⁢                                          ∑                y                            ⁢                                                ∑                  z                                ⁢                                                                  ⁢                                                      H                    ⁡                                          (                                              x                        ,                        y                        ,                        z                        ,                                                  x                          ′                                                ,                                                                              y                            ⁢                                                                                                                                          ′                                                ,                                                  z                          ′                                                                    )                                                        ⁢                                      f                    ⁡                                          (                                              x                        ,                        y                        ,                        z                                            )                                                                                                    +                      n            ⁡                          (                                                x                  ′                                ,                                                      y                    ⁢                                                                                                  ′                                ,                                  z                  ′                                            )                                                          (                  Eq          .                                          ⁢          2.9                )            
Let the RF field (amplitude) g(x′,y′,z′) be measured at a sufficiently large number of points (x′,y′,z′) in a 3D volume space S in the vicinity of the RF source. These points need not be in any regular pattern such as a grid. They need not even be contiguous or close. They could be distributed randomly. But there will be an optimal geometry for the placement of the measurement points which may depend on the target object of study. The minimum number of points is the minimum number of voxels in which the radiation source may be present. But these minimum number of points (x′,y′,z′) must not be a degenerate set, e.g. for certain target objects they should not lie on a surface at constant radial distance. In general, they would have to lie in a 3D volume space that extends substantially along the radial direction. Otherwise the above equation may not be solvable, unless the selected volume V itself has a restricted shape, such as a plane or a line. In the case of a non-degenerate set g(x′,y′,z′), the above equation can be expressed in vector matrix form using conventional notation asg=Hf+n  (Eq. 2.10)
Given measured values of g in a 3D volume space at a set of points (x′,y′,z′), it has been verified through simulation experiments that the above equation can be solved to obtain f. The effect of noise is reduced by using a suitable optimization method in the prior art.